<?php 

namespace App\MoonBundle\Service;

class MoonPhaseService
{
    private $timestamp;
    private $phase;
    private $illum;
    private $age;
    private $dist;
    private $angdia;
    private $sundist;
    private $sunangdia;

    private $synmonth;

    private $quarters = null;

    function __construct( $pdate = null ) {
        if( is_null( $pdate ) )
            $pdate = time();

        /*  Astronomical constants  */
        $epoch = 2444238.5;         // 1980 January 0.0

        /*  Constants defining the Sun's apparent orbit  */
        $elonge = 278.833540;       // Ecliptic longitude of the Sun at epoch 1980.0
        $elongp = 282.596403;       // Ecliptic longitude of the Sun at perigee
        $eccent = 0.016718;         // Eccentricity of Earth's orbit
        $sunsmax = 1.495985e8;      // Semi-major axis of Earth's orbit, km
        $sunangsiz = 0.533128;      // Sun's angular size, degrees, at semi-major axis distance

        /*  Elements of the Moon's orbit, epoch 1980.0  */
        $mmlong = 64.975464;        // Moon's mean longitude at the epoch
        $mmlongp = 349.383063;      // Mean longitude of the perigee at the epoch
        $mlnode = 151.950429;       // Mean longitude of the node at the epoch
        $minc = 5.145396;           // Inclination of the Moon's orbit
        $mecc = 0.054900;           // Eccentricity of the Moon's orbit
        $mangsiz = 0.5181;          // Moon's angular size at distance a from Earth
        $msmax = 384401;            // Semi-major axis of Moon's orbit in km
        $mparallax = 0.9507;        // Parallax at distance a from Earth
        $synmonth = 29.53058868;    // Synodic month (new Moon to new Moon)
        $this->synmonth = $synmonth;
        $lunatbase = 2423436.0;     // Base date for E. W. Brown's numbered series of lunations (1923 January 16)

        /*  Properties of the Earth  */
        // $earthrad = 6378.16;             // Radius of Earth in kilometres
        // $PI = 3.14159265358979323846;    // Assume not near black hole

        $this->timestamp = $pdate;

        // pdate is coming in as a UNIX timstamp, so convert it to Julian
        $pdate =  $pdate / 86400 + 2440587.5;

        /* Calculation of the Sun's position */

        $Day = $pdate - $epoch;                             // Date within epoch
        $N = $this->fixangle((360 / 365.2422) * $Day);      // Mean anomaly of the Sun
        $M = $this->fixangle($N + $elonge - $elongp);       // Convert from perigee co-ordinates to epoch 1980.0
        $Ec = $this->kepler($M, $eccent);                   // Solve equation of Kepler
        $Ec = sqrt((1 + $eccent) / (1 - $eccent)) * tan($Ec / 2);
        $Ec = 2 * rad2deg(atan($Ec));                       // True anomaly
        $Lambdasun = $this->fixangle($Ec + $elongp);        // Sun's geocentric ecliptic longitude

        $F = ((1 + $eccent * cos(deg2rad($Ec))) / (1 - $eccent * $eccent)); // Orbital distance factor
        $SunDist = $sunsmax / $F;                           // Distance to Sun in km
        $SunAng = $F * $sunangsiz;                          // Sun's angular size in degrees

        /* Calculation of the Moon's position */
        $ml = $this->fixangle(13.1763966 * $Day + $mmlong);             // Moon's mean longitude
        $MM = $this->fixangle($ml - 0.1114041 * $Day - $mmlongp);       // Moon's mean anomaly
        $MN = $this->fixangle($mlnode - 0.0529539 * $Day);              // Moon's ascending node mean longitude
        $Ev = 1.2739 * sin(deg2rad(2 * ($ml - $Lambdasun) - $MM));      // Evection
        $Ae = 0.1858 * sin(deg2rad($M));                                // Annual equation
        $A3 = 0.37 * sin(deg2rad($M));                                  // Correction term
        $MmP = $MM + $Ev - $Ae - $A3;                                   // Corrected anomaly
        $mEc = 6.2886 * sin(deg2rad($MmP));                             // Correction for the equation of the centre
        $A4 = 0.214 * sin(deg2rad(2 * $MmP));                           // Another correction term
        $lP = $ml + $Ev + $mEc - $Ae + $A4;                             // Corrected longitude
        $V = 0.6583 * sin(deg2rad(2 * ($lP - $Lambdasun)));             // Variation
        $lPP = $lP + $V;                                                // True longitude
        $NP = $MN - 0.16 * sin(deg2rad($M));                            // Corrected longitude of the node
        $y = sin(deg2rad($lPP - $NP)) * cos(deg2rad($minc));            // Y inclination coordinate
        $x = cos(deg2rad($lPP - $NP));                                  // X inclination coordinate

        $Lambdamoon = rad2deg(atan2($y, $x)) + $NP;                     // Ecliptic longitude
        $BetaM = rad2deg(asin(sin(deg2rad($lPP - $NP)) * sin(deg2rad($minc))));     // Ecliptic latitude

        /* Calculation of the phase of the Moon */
        $MoonAge = $lPP - $Lambdasun;                                   // Age of the Moon in degrees
        $MoonPhase = (1 - cos(deg2rad($MoonAge))) / 2;                  // Phase of the Moon

        // Distance of moon from the centre of the Earth
        $MoonDist = ($msmax * (1 - $mecc * $mecc)) / (1 + $mecc * cos(deg2rad($MmP + $mEc)));

        $MoonDFrac = $MoonDist / $msmax;
        $MoonAng = $mangsiz / $MoonDFrac;                               // Moon's angular diameter
        // $MoonPar = $mparallax / $MoonDFrac;                          // Moon's parallax

        // store results
        $this->phase = $this->fixangle($MoonAge) / 360;                 // Phase (0 to 1)
        $this->illum = $MoonPhase;                                      // Illuminated fraction (0 to 1)
        $this->age = $synmonth * $this->phase;                          // Age of moon (days)
        $this->dist = $MoonDist;                                        // Distance (kilometres)
        $this->angdia = $MoonAng;                                       // Angular diameter (degrees)
        $this->sundist = $SunDist;                                      // Distance to Sun (kilometres)
        $this->sunangdia = $SunAng;                                     // Sun's angular diameter (degrees)
    }

    private function fixangle($a) {
        return ( $a - 360 * floor($a / 360) );
    }

    //  KEPLER  --   Solve the equation of Kepler.
    private function kepler($m, $ecc) {
        //double e, delta;
        $epsilon = pow(1, -6);
        $e = $m = deg2rad($m);
        do {
            $delta = $e - $ecc * sin($e) - $m;
            $e -= $delta / ( 1 - $ecc * cos($e) );
        }
        while ( abs($delta) > $epsilon );
        return $e;
    }

    /*  Calculates  time  of  the mean new Moon for a given
        base date.  This argument K to this function is the
        precomputed synodic month index, given by:
            K = (year - 1900) * 12.3685
        where year is expressed as a year and fractional year.
    */
    private function meanphase($sdate, $k){
        // Time in Julian centuries from 1900 January 0.5
        $t = ( $sdate - 2415020.0 ) / 36525;
        $t2 = $t * $t;
        $t3 = $t2 * $t;

        $nt1 = 2415020.75933 + $this->synmonth * $k
                + 0.0001178 * $t2
                - 0.000000155 * $t3
                + 0.00033 * sin( deg2rad( 166.56 + 132.87 * $t - 0.009173 * $t2 ) );

        return $nt1;
    }

    /*  Given a K value used to determine the mean phase of
        the new moon, and a phase selector (0.0, 0.25, 0.5,
        0.75), obtain the true, corrected phase time.
    */
    private function truephase($k, $phase){
        $apcor = false;

        $k += $phase;               // Add phase to new moon time
        $t = $k / 1236.85;          // Time in Julian centuries from 1900 January 0.5
        $t2 = $t * $t;              // Square for frequent use
        $t3 = $t2 * $t;             // Cube for frequent use
        $pt = 2415020.75933         // Mean time of phase
             + $this->synmonth * $k
             + 0.0001178 * $t2
             - 0.000000155 * $t3
             + 0.00033 * sin( deg2rad( 166.56 + 132.87 * $t - 0.009173 * $t2 ) );

        $m = 359.2242 + 29.10535608 * $k - 0.0000333 * $t2 - 0.00000347 * $t3;          // Sun's mean anomaly
        $mprime = 306.0253 + 385.81691806 * $k + 0.0107306 * $t2 + 0.00001236 * $t3;    // Moon's mean anomaly
        $f = 21.2964 + 390.67050646 * $k - 0.0016528 * $t2 - 0.00000239 * $t3;          // Moon's argument of latitude
        if ( $phase < 0.01 || abs( $phase - 0.5 ) < 0.01 ) {
           // Corrections for New and Full Moon
            $pt +=  (0.1734 - 0.000393 * $t) * sin( deg2rad( $m ) )
                    + 0.0021 * sin( deg2rad( 2 * $m ) )
                    - 0.4068 * sin( deg2rad( $mprime ) )
                    + 0.0161 * sin( deg2rad( 2 * $mprime) )
                    - 0.0004 * sin( deg2rad( 3 * $mprime ) )
                    + 0.0104 * sin( deg2rad( 2 * $f ) )
                    - 0.0051 * sin( deg2rad( $m + $mprime ) )
                    - 0.0074 * sin( deg2rad( $m - $mprime ) )
                    + 0.0004 * sin( deg2rad( 2 * $f + $m ) )
                    - 0.0004 * sin( deg2rad( 2 * $f - $m ) )
                    - 0.0006 * sin( deg2rad( 2 * $f + $mprime ) )
                    + 0.0010 * sin( deg2rad( 2 * $f - $mprime ) )
                    + 0.0005 * sin( deg2rad( $m + 2 * $mprime ) );
            $apcor = true;
        } else if ( abs( $phase - 0.25 ) < 0.01 || abs( $phase - 0.75 ) < 0.01 ) {
            $pt +=  (0.1721 - 0.0004 * $t) * sin( deg2rad( $m ) )
                    + 0.0021 * sin( deg2rad( 2 * $m ) )
                    - 0.6280 * sin( deg2rad( $mprime ) )
                    + 0.0089 * sin( deg2rad( 2 * $mprime) )
                    - 0.0004 * sin( deg2rad( 3 * $mprime ) )
                    + 0.0079 * sin( deg2rad( 2 * $f ) )
                    - 0.0119 * sin( deg2rad( $m + $mprime ) )
                    - 0.0047 * sin( deg2rad ( $m - $mprime ) )
                    + 0.0003 * sin( deg2rad( 2 * $f + $m ) )
                    - 0.0004 * sin( deg2rad( 2 * $f - $m ) )
                    - 0.0006 * sin( deg2rad( 2 * $f + $mprime ) )
                    + 0.0021 * sin( deg2rad( 2 * $f - $mprime ) )
                    + 0.0003 * sin( deg2rad( $m + 2 * $mprime ) )
                    + 0.0004 * sin( deg2rad( $m - 2 * $mprime ) )
                    - 0.0003 * sin( deg2rad( 2 * $m + $mprime ) );
        if ( $phase < 0.5 )     // First quarter correction
            $pt += 0.0028 - 0.0004 * cos( deg2rad( $m ) ) + 0.0003 * cos( deg2rad( $mprime ) );
        else    // Last quarter correction
            $pt += -0.0028 + 0.0004 * cos( deg2rad( $m ) ) - 0.0003 * cos( deg2rad( $mprime ) );
            $apcor = true;
        }
        if (!$apcor)    // function was called with an invalid phase selector
            return false;

        return $pt;
    }

    /*  Find time of phases of the moon which surround the current date.
        Five phases are found, starting and
        ending with the new moons which bound the  current lunation.
    */
    private function phasehunt() {
        $sdate = $this->utctojulian( $this->timestamp );
        $adate = $sdate - 45;
        $ats = $this->timestamp - 86400 * 45;
        $yy = (int) gmdate( 'Y', $ats );
        $mm = (int) gmdate( 'n', $ats );

        $k1 = floor( ( $yy + ( ( $mm - 1 ) * ( 1 / 12 ) ) - 1900 ) * 12.3685 );
        $adate = $nt1 = $this->meanphase( $adate, $k1 );

        while (true) {
            $adate += $this->synmonth;
            $k2 = $k1 + 1;
            $nt2 = $this->meanphase( $adate, $k2 );
            // if nt2 is close to sdate, then mean phase isn't good enough, we have to be more accurate
            if( abs( $nt2 - $sdate ) < 0.5 )
                $nt2 = $this->truephase( $k2, 0.0 );
            if ( $nt1 <= $sdate && $nt2 > $sdate )
                break;
            $nt1 = $nt2;
            $k1 = $k2;
        }

        // results in Julian dates
        $data = array(
            $this->truephase( $k1, 0.0 ),
            $this->truephase( $k1, 0.25 ),
            $this->truephase( $k1, 0.5 ),
            $this->truephase( $k1, 0.75 ),
            $this->truephase( $k2, 0.0 )
        );

        $this->quarters = array();
        foreach( $data as $v )
            $this->quarters[] = ( $v - 2440587.5 ) * 86400; // convert to UNIX time
    }

    /*  Convert UNIX timestamp to astronomical Julian time (i.e. Julian date plus day fraction).  */
    private function utctojulian( $ts ) {
        return $ts / 86400 + 2440587.5;
    }

    private function get_phase( $n ) {
        if( is_null( $this->quarters ) )
            $this->phasehunt();

        return $this->quarters[$n];
    }

    /* Public functions for accessing results */

    function phase(){
        return $this->phase;
    }

    function illumination(){
        return $this->illum;
    }

    function age(){
        return $this->age;
    }

    function distance(){
        return $this->dist;
    }

    function diameter(){
        return $this->angdia;
    }

    function sundistance(){
        return $this->sundist;
    }

    function sundiameter(){
        return $this->sunangdia;
    }

    function new_moon(){
        return $this->get_phase( 0 );
    }

    function first_quarter(){
        return $this->get_phase( 1 );
    }

    function full_moon(){
        return $this->get_phase( 2 );
    }

    function last_quarter(){
        return $this->get_phase( 3 );
    }

    function next_new_moon(){
        return $this->get_phase( 4 );
    }
}